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| \ | INTRODUCTION 
The Naval Ocean Systems Center (NOSC) was tasked by NAVELEX 612 
under the subject of empirical detection analysis to develop techniques 
= for estimating detection probabilities. Computer Sciences Corporation 
>-tescy was subsequently tasked by NOSC to develop data bases and nonlinear 
regression analysis techniques to describe the acoustic detection process 
as a function of signal excess. 


Statistical techniques, primitive empirical data bases and ancillary 
computer programs used to relate the empirical data to models were de- 
veloped previously at NOSC during other analyses of acoustic detection 
probabilities. These resources and the UCLA Department of Biomathematics 

C {BMBPY Statistical Software were made available as a baseline for the 
nonlinear regression analysis methodology used in this investigation. 
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PROBLEM STATEMENT 


Technical Requirements 


The general problem addressed in this task is to establish a meth- 
odology for the determination and use of parameters that describe acoustic 
detection as a function of signal excess. The statistical technique to 
be used is nonlinear regression analysis using the computer programs 
available inthe UCLA BMDP statistical programs library available on the 
UNIVAC 1100/82 system at NOSC. 


Specific subtasks assigned and addressed in this report are to: 


1. Establish the viability of the current version of the BMDP 
o nonlinear regression software. 


2. Prepare data bases for input to the nonlinear regression 
analysis programs. 


3. Write computer programs to implement the Gaussian distribution 
function in the nonlinear model. 


4. Perform nonlinear regression analysis with supplied data sets. 


Approach s 
The derivation of the regression model used in this study is based ~ 
on the sonar equation and the concept of random acoustic fluctuations. j 


The passive sonar equation may be stated in decibel form as: 


a 
SE = SL - TL - AN + AG - RD 
where 
SE is signal excess . 


SL is source level 
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TL is transmission loss 
AN is ambient noise level , 
wn { 
wh AG is array gain 


RD is recognition differential (detection threshold) 
The term performance index (PI) is defined as follows: 
PI = TL + (AN - AG) 


Substituting PI into the sonar equation gives 


Mm. 


SE = (SL - RD) - PI 


It is known that there are random acoustic fluctuations which cause 
. the signal excess to vary over time even though the sonar equation param~ 
| eters are held constant and that these fluctuations have a standard 
: deviation about the mean of approximately 8 dB. The sonar equation can 
be rewritten in a form normalized to the standard deviation of signal 
excess (a). 


Linear Regression Model 


Considering that PI is the independent variable which varies as a 
‘ target moves about within tracking range of an array and that the proba- 
ce bility of detecting and holding a target depends on the signal excess, 
; the above form of the sonar equation suggests a simple linear slope- 
hae intercept regression model. 


Y=P, +P X+e 
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where 


Y is the inverse cumulative Gaussian of (l-fractional holding time) 


X is the performance index 
Py = (SL - RD)/o is the Y-intercept of the regression line 
Po = 1/o is the slope of the regression line 


¢ is the error term 


The linear regression model may be extended to include empirical 
data from several different targets by the use of dummy variables. Then 
the equation becomes: 


Y =P, + PoX, + P3k, 2? Xgt ... + PX, t+ 


1 eo 4 


where 


Ko ses Xn are coded as 1 to indicate the presence of specific 
additional targets data and 0 otherwise. 


This arrangement as a multiple linear regression model provides for 
separate parallel regression lines with different intercepts for each 
target contributing data to the sample but a common slope for all re- 
gression lines. The advantage of multiple regression analysis is a 
better estimate of the slope due to increased sample size which for 
practical cases greatly overrides the disadvantage of giving up an 
additional degree of freedom for each dummy variable added. 


Nonlinear Regression Model 


Previous regression analyses performed at NOSC used the linear 
regression model described in the previous paragraphs. However, the 
current task is to perform the regression analyses using nonlinear 
regression, that is, a model of the form: 


Se Ce ee Re Sy Bs MRT eS ngs Ee Me ee ee ee a EN Ne See ieee Ke ee gee: 


Ot. 


| Y=s=1-F(Z) +e 


where 


1 Z -(t2)/2 
F(Z) = e dt; the cumulative Gaussian distribution 


bel i 
es + function 
Z = Pr + Pox + P 4X5 + scat Panel 
Y is the probability of detecting/holding a contact (assumed to be 
a nonlinear function of Pi» X;) 
ad ¢ is the error term 
z This model will permit the use of a different weighting function 
; than that used with the linear model which effectively excluded some 
| data points by zero weighting. Hopefully, a more accurate determination 
of regression parameters will result. The BMDP Statistical Software 
(references 1 and 2) available at NOSC contains computer programs suit- 
able for implementing this nonlinear regression model. 
a 
Weighting Function 
There is no basis to assume that the variance of the error term 
= in the regression equation is homogeneous, therefore it is necessary to 
2 use case weighting (reference 1, p 302) inversely proportional to the 
variance. The weighting function provided for this study is: 
oe BF 
W. = 
ee where 
a , 
W, is the weighting assigned to a datum 
we 5 
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T, is the available holding time for the datum 


F(Z.) is the cumulative Gaussian distribution function evaluated for 
the datum 


The weighting function is used in the nonlinar regression algorithm 
as part of the least squares criterion for a best fit to the observed 
data. That is, when determining the parameters for regression line equa- 
tions, the term 


N 

a 2 
» Wi(¥5 = YF) 
i=l 


is minimized where: 


i is the case index 
W is the weight for a particular case 
Y; is the observed fractional holding time 


VF is the predicted fractional holding time 
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RESULTS AND DISCUSSION 


Nonlinear Regression Computer Program Viability 


BMOP Software, 1977 Version 


At the commencement of this study, the 1977 version of the UCLA 
BMDP Statistical Software library (reference 2) was operational on the 
UNIVAC 1100/82 at NOSC and the 1981 version of the programs (reference 1) 
was undergoing installation. Initial program evaluation was accomplished 
using the 1977 versions of both the P3R and the PAR nonlinear regression 
analysis programs (reference 2, p 464). First the example problems pre- 
sented with the program description in the BMDP documentation were 
executed and the outputs verified with the published outputs. The examples 
were found to execute as advertised in the BMDP document. Next an example 
was taken from a textbook (reference 3, example 10.4) which had been 
worked out using a nonlinear regression algorithm different than the 
one implemented in the BMDP programs. The nonlinear function implement- 
ing this example was programmed and collected with both the P3R and PAR 
programs. The regression equation parameters obtained by executing P3R 
and PAR matched those in the textbook example. Finally, three sets of 
data were taken from prior linear regression analyses and processed 
through the P3R and PAR software programmed for the cumulative Gaussian 
nonlinear regression model with iteratively reweighted least squares 
described in previous paragraphs. Results from executing the programs 
yielded parameters which compared favorably with the linear regression. 
Appendix A lists the formatted outputs from the 1977 version of the non- 
linear regression analysis programs. 


Comparison of P3R and PAR Programs 


The BMDP P3R and PAR programs both perform nonlinear regression by 
the minimizing of least squares criterion. The algorithms differ opera- 
tionally in that P3R requires the first partial derivatives of the nonlinear 


function with respect to each parameter in the regression equation while 
PAR does not require them. Parameter estimates output from the two pro- 
grams were generally in good agreement for the cases tested but PAR re- 
quired more iterations to converge to a solution than P3R. Since the 
partial derivatives and satisfactory initial estimates can be specified 
for the cumulative Gaussian function, the P3R program was selected as 
the best algorithm to use for this study. 
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rn BMDP Software, 1981 Version 


When the 1981 version of the BMDP P3R program was installed on the 
UNIVAC 1100/82 the same checkout procedures as executed for the 1977 
version were used. The results from the two versions were found to be 
essentially in agreement. The 1981 version of the program contains some 
more advanced features than the 1977 version and appears to converge with 
fewer iterations so the 1981 version of the P3R program is the best pro- 
gram to use for the nonlinear regression analyses assigned to this task. 
Appendix B lists the outputs from the execution of the 1981 version of the 
P3R nonlinear regression analysis program. 


Evaluation of Cumulative Gaussian Distribution Function “ 


Initial efforts with the 1977 BMDP programs used the single precision 
subroutine MDNOR available in the IMSL mathematical subroutine library on 
the UNIVAC 1100/82 at NOSC (reference 4). The BMDP P3R and PAR programs 
both operate in double precision but the IMSL double precision subroutine 
MDNORD is not implemented in the IMSL-8 library at NOSC. It became appar- 
ent during runs with the nonlinear regression model that a double preci- 
sion subroutine for evaluation of the cumulative Gaussian distribution = 
function would be required because the dynamic range of the input argu- 
ment during the regression program iterations sometimes exceeded the 
capability of the single precision subroutine to calculate the function. 


MONORD Subroutine Algorithm 


A double precision subroutine, called MDNORD after the nonimplemented 
IMSL subroutine, was developed for use with the nonlinear regression analy- 
= sis model. This subroutine is based on the recursive evaluation of continuing 
fraction expressions which estimate the area under the normal or Gaussian 
distribution function (reference 5, arts. 26.2.14 and 26.2.15). The 
algorithm developed for subroutine MDNORD uses the art. 26.2.14 equation 
for arguments where the magnitude of x is greater than three, i.e., 


Ax) = 2x4 ESSE dS |x > 3 


o 


and when the argument is equal to or less than three the art. 26.2.15 
equation is used. 


a 2 ce mee 
| ax) =e 20} PERM be el cs 
where 
| x is the abscissa of the Gaussian distribution function expressed as 


the number of standard deviations from the mean. 
Q(x) is the complement of the area under the Gaussian (1l-area). 


Z(x) is the ordinate of the Gaussian distribution function. 


This dual algorithm was used because testing showed that the number 
: of terms required for the art. 26.2.14 equation to reach double precision 
oa accuracy increased as the input argument x became smaller whereas the 
art. 26.2.15 equation exhibited the reverse effect with the number of 
terms required increasing as x became larger. A convenient branching 
point was reached at a value of x = 3 standard deviations from the mean. 


The FORTRAN program listing and test outputs for subroutine MDNORD 
are included in Appendix C of this report. The source and relocatable 
code are located in PASS*NRL.MDNORD on the UNIVAC 1100/82. 


Continuing Fraction Evaluation 


A convenient theorem on the evaluation of continued fractions (refer- 


ence 5, art. 3.10) was implemented in subroutine MDNORD. This theorem 
provides a matrix multiplication technique for evaluation of the nth 
tional of the continued fraction. 


A a, a, a a 
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frace 


A, and B, are by definition terms expressing the numerator and denominator 


of the nth fractional which are used in stating the matrix form of the 
theorem. 
A Anel An-2 DF 


n 


By Bal By-2 ay 


Calculation of the nth fractional can be done recursively once the 
two matrices on the right side of the equation are initialized. The 
MONORD program uses a recursive calculation loop for t. = (A,)/8, which 
ae until the difference between consecutive terms is less than 
1% 10° 7%; 


P3RFUN and FUN Subroutines 


The Cumulative Gaussian Nonlinear Regression function described 
previously is implemented as subroutine P3RFUN. This subroutine is 
written in FORTRAN and is listed in Appendix C. The logic is set up to 
evaluate the expression 
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where 


Z = Py + Pax, + PX, 5 ee i Pixel 


and the first partial derivatives of Y with respect to the nonlinear 
parameters Py> Pos wee, P. (see reference 6, art. 67, differentiation 
under the integral) 


acs 


dY(Z dY(Z 
ae 
2...n-1 2...N d 1 


as well as the weighting function 


2 
Ty 


Wi * TD FUN) F(Z.) 
i * T= FJ) FU, 


each time the BMDP P3R nonlinear regression program calls P3RFUN. 


Each iteration of the nonlinear regression analysis algorithm causes 
a call to P3RFUN for each case in the input data set. Thus the total 
number of calls to P3RFUN during a nonlinear regression analysis run 
is the number of iterations times the number of cases. Program branches 
are provided to write a debug aid file on logical unit 20 if the BMDP 
P3R input NUMBER jis set to a value of 3. The data written in the debug 
file is the sequential number of the call to P3RFUN, the case number, the 
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3 1 21 3°2 
os distribution curve and the value of the weighting function. This feature 


a can create a large amount of data and should be used only under controlled 
conditions to observe the variations in the arguments when a close look 

is needed at the dynamics of the nonlinear regression algorithm. Normally, 
the input NUMBER is set to a value of 2 for iteratively recalculated weight- 
ing. If desired, any constant weighting input is used via a branch 

which activates when NUMBER is set to 1. 


ra value of Z = P, + PX, + P3X, +... + PAX 1, the area under the Gaussian { 


ise 
Paes ee 


Very large values of the weighting function sometimes occur for 
extreme values of the input argument Z to the P3RFUN program. These can 
be seen to cause the factor Y(Z)[1 - Y(Z)] to approach zero for either 
wg Y(Z) approaching 1 or Y(Z) approaching 0 which corresponds to a Z value 
: near one or the other tails of the Gaussian distribution function. When- 
ever the factor [Y(Z)][1 - Y(Z)] becomes zero by exceeding precision 
limitations the factor is reset to 1.0 x 10738 in order to avoid division . 
by zero. In this event, the sequential call number, case number, Z value, - 
area under the Gaussian and the factor are all written to a file on 
logical unit 21. In addition, the condition word is set so that the 
executive can sense the condition and the file can be printed out. 


a Subroutine FUN is the equivalent of subroutine P3RFUN without the 
evaluation of the derivatives. This subroutine is used with the PAR 

3 program in the same way that P3RFUN is used with the P3R program. The 
"\ FORTRAN source code for both P3RFUN and FUN are included in Appendix C 
of this report. 


Collection of Executable Programs 


The MAP directives used to collect the cumulative Gaussian nonlinear 
regression model programs are listed in Appendix C. The executable pro- 
grams are located in files PASS*NLR.P3R81, PASS*NLR.P3R77 and ie 
PASS*NLR.PAR77. The BMDP programs provided in the NOSC implementation 3 
on the UNIVAC 1100/82 have computer programs available to assist the user 
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in constructing the MAP as well as the externally provided nonlinear func- 
tion. These programs are invoked by using the appropriate EXEC 8 command 
from the following list: 


@N*BMDP77 . 3RBUILDF 
@N*BMOP 77. ARBUILDF 
@N*BMOP81.BUILDFUN/P3R 
@N*BMDP81.BUILDFUN/PAR 


The user then inputs the external function FORTRAN code excluding the 
DIMENSION or VARIABLE typing and the RETURN,END. 


Output from the above processors is written to logical unit 8 so 
the program may be collected (@MAPed) and executed in file TPF$. by 
simply using the GADD 8 command. However, if repeated executions of the 
absolute element are desired, as in this project, it is wasteful of time 
and computer resources to collect the program each run. Therefore the 
MAP commands were extracted from the logical unit 8 file by means of the 
MED text editor and saved for use as a separate set of directives. 

The user may also save the FORTRAN code for the external function in a 
similar way if desired. 


Input Data Base Construction 


Empirical detection data was taken from two sources and reformatted 
for input to the BMDP nonlinar regression analysis programs. The input 
data base for the PAR and P3R programs is constructed as symbolic elements 
in the FILE.ELEMENT format of the UNIVAC 1100/82 EXEC 8. Two programs 
were written to reformat data files from either the linear program input 
format or the primitive data format. 


FILEPROC Program 


The FILEPROC program is designed to accept data in the format used 
for the linear regression analysis program and to reformat it for input 
to the nonlinear regression analysis programs. Operation of the program 
is interactive in demand mode at a terminal. The EXEC 8 command 
@XQT PASS*NLR.FILEPROC runs the program. The user should have input data 
available in either files or elements of files so that it can be added 
to the runstream at the terminal by the @ADD command when solicited. The 
output is written to a formatted data file on logical unit 8 and the data 
can be moved about from the unit 8 file by the MED processors. The FORTRAN 
code and an example execution of the FILEPROC program are included in 
Appendix D. 


FILPROCS Program 


The FILPROCS program is similar to the FILEPROC program in that it 
reformats data for input to the nonlinear regression analysis programs. 
However, it receives data from a number of different input files (one 
for each sound projector or source recorded at an array) and combines the 
data into a single file. Therefore, the user must have at hand the set 
of files, previously extracted from primitive data, which represents the 
whole data set to be used for a nonlinear regression analysis run. The 
program is executed by the command @XQT PASS*NLR.FILPROCS and the file 
inputs are solicited interactively. It should be noted that the runstream 
entry @EOF is necessary after each data file is added as well as after 
all data is completed. The @EOF may be inserted in the input data file 
as shown in the FILPROCS example of Appendix D. These @EOFs are not a 
part of the files after they are extracted from the primitive data, j.e., 
they were inserted by means of the MED processor for convenience. Output 
from the FILPROCS program is written to a data file on logical unit 8. 
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Post Processing Program - SXPROC 


The BMDP P3R program data save files do not provide all quantities 
needed for the analyses of this project. Therefore a post processing 
program was developed to read selected data from saved printout files 
and process these data for display, particularly display on the basis 
of the signal excess. At this writing, SXPROC performs two functions, 
namely, collection of holding time data on the basis of signal excess and 
calculation of the parameter covariance matrix. 


Signal Excess Versus FHT 


The interpretation of the regression equation parameters as related 
to the sonar equation is as follows: 


a= 1/P, is the standard deviation of the signal excess about 
the mean. 


SL - RD =0 Py is the source level - recognition differential 
difference. 


SX = (SL - RD) - PI is the signal excess. 


where 


PI is the performance index which is the independent variable in the 
regression analysis. 


The holding time and available holding time are transformed from 
the 1.5 dB PI domain intervals and accumulated in the corresponding SX 
domain 1.5 dB intervals. After the data is accumulated for each source 
in this manner, it is summed over al] sources to form the holding time 
data for the sensor. Then the nonholding time can be formed and com- 
bined with a count of the gains (number of source holding periods follow- 
ing a previous source loss) for each SX interval and used in likelihood 
ratio analyses of the data. 

15 
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Parameter Covariance Matrix 


The covariance matrix of the regression equation parameters is useful 
for analyzing the confidence intervals for the prediction equations. It 
is calculated from the BMDP P3R estimates of the standard deviations of 
the parameters and the correlation matrix for all the parameters using 
the relationship: 


cov(P., P,) = S, S; Mig 


S; is the standard deviation estimated for a parameter P. 
and 


r.. is the correlation coefficient between two parameters P. and P 


iJ a” 


The covariance matrix is stored and processed in lower left triangular 
format and a formatted printout of the matrix is available. 


SXPROC Operation 


The SXPROC program is designed for the FORTRAN NAMELIST type of 
input. The inputs available currently are: 


COVAR - A logical variable set to T or F to initiate the covariance 
matrix calculation and display. 


SIGEX - A logical variable set to T or F to initiate the FHT versus 
signal excess processing and display. 


FILNAM - A character variable for inputting the data file name (up 
to 12 characters). 


ELEM - A character variable for inputting an element name (up to 
12 characters). 
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VERS - A character variable for inputting a version name (up to 12 
characters). 


UNIT - The logical unit to use for readin of the data file (default 
is LU 8). 


STOP - A logical variable set to T for program stop, otherwise another 
set of namelist data is expected. 


The execution of SXPROC is by the EXEC 8 command @XQT PASS*NLR.SXPROC. 


The method of constructing a file for input is to breakpoint the 
print files for a P3R nonlinear regression run as illustrated in the 
example runstream for P3R81 in Appendix B. The MED processor is then 
used to save the printout file as an element/version of the intermediate 
file needed for input to SXPROC. Typically, element names might be 
selected to identify arrays and versions to identify variations of con- 
ditions at that array such as season. An example runstream for SXPROC 
using elements built from the P3R81 example is shown in Appendix E along 
with the formatted printouts resulting from the run. The FORTRAN source 
code for the SXPROC programs is also included in Appendix E. 


Post Processing Program - NLRPLT 


The NLRPLT program is designed to plot the nonlinear regression 
analysis prediction curves. It uses the saved printout files as input 
just as the SXPROC program does. Two different kinds of plots are 
plotted for each regression analysis run. The first kind of plot dis- 
plays the prediction for each source contributing data to the analysis 
at a particular array with all curves on a single plot. The observed 
data, from which the predictions result, are shown in the same graph 
marked with X's of varying size. The size of the X's marking data are 
proportional to the square root of the case weights. Therefore it can 
be seen which points have greater or lesser amounts of influence on the 
regression curve. The second kind of plot displays the prediction 


17 


Se ee ee ae, et ge et de re i ee To com Pee aA 
fig Seg gee Tg et he ar e at tala : 3 : F : 


curve for one particular source along with the curves bounding the 95 
percent confidence region for the predicted curve. A separate plot is 
made for each source in the data set. 


The NLRPLT program is also capable of printing out the source level- af 
ee recognition differential differences for each source and the covariance 
matrix for the parameters of the regression equation. 
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NLRPLT Operation 


: The NLRPLT program is designed for the FORTRAN NAMELIST type of 
= input similar to the SXPROC program. The inputs available are: 


{ I " 2 


FILNAM - A character variable for inputting the data file name 
(up to 12 characters). 


ELEM - A character variable for inputting the element name in the : 
data file (up to 12 characters). n 


VERS - A character variable for inputting a version name (up to 12 
characters). 


UNIT - The logical unit to use for readin of the data file (default cs 
is LU 8). 


COVAR - A logical variable set to T or F to indicate display of 
the covariance matrix. 


: LEVELS - A logical variable set to T or F to indicate display of 
- source level-recognition differential display. 


PLOTS - A logical variable set to T or F to indicate plot displays. 


STOP - A logical variable set to T or F for program stop, otherwise 
another set of namelist data is expected. 


The execution of the NLRPLT program is by the EXEC 8 command 


@XQT PASS*NLRPLT.NLRPLT. 
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Input file construction for NLRPLT is identical to that described 
for the SXPROC program (refer to previous paragraphs on SXPROC). An 
example runstream with the resulting plot and print outputs are included 
as Appendix F of this report. The FORTRAN source code and the MAP collec- 
tion directives are also included. 


Maximum Likelihood Criterion 


Program P3R may be used with a maximum likelihood calculation to fit 
the nonlinear regression function to the data rather than the least 
squares criterion mentioned under the weighting function description. 

The method recommended (see reference 1, p 315) for recursive recalcula- 
tion of the weighting function when fitting a curve such as the Gaussian 
to a data set is as follows: 


Turn off the convergence criterion by using the following 
inputs to P3R. 


CONVERGENCE IS -1.0. 
HALVING IS 0. 


Rescale the standard deviations of the parameters to convert 
them to standard errors by the input. 


MEANSQUARE IS 1.0. 


This method used with a maximum of 15 iterations gives the best 
results for nonlinear regression parameters. A sample runstream and P3R 
program outputs are included in Appendix B as run P3R81A. 
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CONCLUSIONS AND RECOMMENDATIONS 
Conclusions 


Based on the results of computer runs during this task it is 
concluded that: 


e The P3R program is the better nonlinear regression analysis 
program for fitting the cumulative Gaussian function to ob- 
served data for this task. 


e The 1981 version of the P3R program gives results essentially 
the same as the 1977 version. 


e Data observations indicating the absence of detection as well 
as data indicating detections are effectively included in the 
nonlinear regression. 

Recommendations 
It is recommended that: 
e Future nonlinear regression analysis runs use the 1981 version 


of the UCLA BMDP P3R program and the maximum likelihood 
criterion. 
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